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Abstract. 

We compute numerical simulations of spherical collapse triggered by a slow increase in external pressure. We 
compare isothermal models to models including cooling with a simple but self-consistent treatment of the coupling 
between gas, grains and radiation field temperatures. The hydrostatic equilibrium appears to hold past the 
marginally stable state, until the collapse proceeds. The last hydrostatic state before collapse has a lower central 
gas temperature in the centre due to the enhanced coupling between gas, grains and radiation field. This results 
in slightly lower pressure gradients in the bulk of the envelope which is hence slightly more extended than in the 
isothermal case. Due to the sensitivity of the collapse on these initial conditions, protostellar infall velocities in 
the envelope turn out to be much slower in the case with cooling. 

Our models also compute the radiative transfer and a rather large chemical network coupled to gas dynamics. 
However, we note that the steady-state chemisorption of CO is sufficient to provide an accurate cooling function 
of the gas. This justifies the use of post-processing techniques to account for the abundance of observed molecules. 
Existing observations of infall signatures put very stringent constraints on the kinematics and temperature profile 
of the class protostar IRAM 04191-1-1522. We show that isothermal models fail to account for the innermost 
slow infall motions observed, even with the most hydrostatic initial conditions. In contrast, models with cooling 
reproduce the general shape of the temperature profile inferred from observations and are in much better agreement 
with the infall signatures in the inner 3000 AU. 
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1. Introduction Il999|) have paid extreme attention to the energy budget, 

r- 1 I -I using accurate models of radiative transfer coupled with 

We revisit here the n ume rical work or ILarsonl II1969I I. riji,-i j. 

r— -1 I 1 r— o T 1 I J r£ complicated chemical networks designed to account 

IMasunaea et al.l II1998I) and lMasunaea fc Inutsukal II2000I1 r .> i i f .i • r j. r .i 

I —ta— -1 1 i4/ L -la- ta 1 la/ for the abundances or the mam cooling agents oi the 

to probe the influence of cooling and chemistry on the col- • 4. j. n j- tt • • 1 4-1, 1 j 1 

T ^ n 1 interstellar medium. Using more simple thermal models, 

lapse ol protostellar and prestellar dense cores. We com- UTTrZ 71 /^oooo^ i, 1 • i- ^ j ^, j i j.- 

, , r ^1 ICralli et al.l H2002ii have also investigated hydrostatic 

pare our results to observations of a Class protostar. I.e., . • n 4. ui 4. 4. 

„ „ , „ equilibria near the marginally stable state, 

a collapsing protostellar core so young that i t still retains 

J 4. , • , r -4. • •4.- 1 J-4-- n T~' I n In the present wo r k, we use a similar approach to 

a det ailed picture 01 its initial conditions IIAndre et alJ 1— f— — 1 j ' 1 , ^ , 

Ij^gggj^ IMasimaga fc Inutsukal H2000I) . We use a simplihed treat- 

T ,1 1,1- r 1 • 1 11 1 r ment of the radiation transfer, which is found to give very 
Isothermal studies of sph erical collapse have to- ... m, • „ , ,1 ^t^tt • 
, ,„ . , Jni l^^f^r,-, T ir>^^ similar results, ims allows us to spend the CFU time on 
cused on self- similar solutions llshul 11977: La rson T969t , ., , , . , , , ttt 
^"^^^ionK • •,• 1 a more detailed chemical network. We are hence able to 



iBouauet et al.l Il985t iBlottiau et all Il988l ). initial con , , , , .,, , r 

Jt^ i D r^i 1- J I1 ^ — j2 1 1 modcl thc molecular cooling with a degree oi accuracy 

ditions H l'oster & Chevalier! 119931) and magnetic faeids . i i i • i i a i 

nr^ — ^ ■ ,„„. T- 1f^f,o^ closc to the hydrostatic models. Also, we use a moving 

llBasu fc M ouschovias 1994: Li 1998). ., , , ' . , , 

' • 1 • r 1 1 • grid algorithm which describes the accretion shock with a 

In contras t, simul atio ns ot hydrostatic clouds i i 4_4_ i 4_- 

rr — ; — 71 r j 177- , „ '' — 1 ITTTTl much better resolution. 

llDe Jong et alJ Il980t iBoland fc De J ong 1984"; 

IChieze fc Pineau des ForetA Il987t iNeiad fc Wagcnblast ^he aim of this paper is to show how the use of a 

more accurate model for the gas and dust microphysics 

Send offprint requests to: lesaffre@ast.cam.ac.uk can influence the dynamics of the protostellar collapse. We 
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briefly present isothermal models (Section 2) before com- 
paring them to the results of computations including radi- 
ation transfer and using a full chemical network (Section 
3). We emphasise the importance of the coupling model 
between gas and radiation through the grains which al- 
lows us to compute a realistic profile for the temperature. 
This profile is responsible for slightly lower pressure gradi- 
ents which make the collapse milder than in the isothermal 
case. We investigate the influence of a few parameters and 
in particular show that the out-of-equilibrium chemistry 
has little influence on the cooling of the gas. We finally 
compare our models to observations of the young low mass 
class protostar IRAM 04191 (Section 4). We find that 
isothermal models fail to reproduce the slow observed in- 
fall velocities. On the other hand, models with cooling are 
within the observational uncertainties inside a radius of 
around 3000 AU. We discuss and sum up our results in 
sections 5 and 6. 



2. Isothermal models 

2.1. Initial conditions 

We start our collapse calculations from a stable hydro- 
static state. We then trigger the gravitational instability 
by a slow exponential increase of the external pressure 
Pext- When the critical pressure is reached, the collapse 
evolves on a few free-fall time scales, until the birth of the 
protostar. A jump in the mass-to-radius ratio at the cen- 
tre characterises this particula r instant, wh ich wc take as 
the origin of times t = (as in lFoster fc Che valier 1993) . 
In the present work, we refer to times t < as the prestel- 
lar phase and t > as the protostellar phase, in agree- 
ment with the terminology used by most observers. Three 
parameters specify a simulation : the total mass Af, the 
temperature T, and the time scale for the pressure in- 
crease tp =dt/dln(poxt)- We use = 2.33 amu for the 
mean molecular weight . 

Be cause the eq uations can be put in an adimensional 
form l)Foster fc Chevalier 1993') only the parameter tp is 
relevant. We chose to fix M = 1.7 Mq and T = 10 K 
as good estimates for the object we plan to compare to 
(see section ^21- The free-fall time-scale of the marginal 
Bonnor-Ebert sphere with these parameters is iff = 1.3 x 
10^ yr. 



2.2. Method 

The use of a Lagrangian mesh allows us to easily treat 
the huge dynamical range in densities and scales in- 
volved in the collapse without introducing advection er- 
rors. However, it is not able to evolve through the cen- 
tral singularity. We therefore use a moving grid algorithm 
that keeps the mesh on a logarithmic spacing in radius 
with a fixed radius for the central zone when it gets below 
10~* pc. To check this method, we compared our prestellar 
results with the output of a Lagrangian code. The com- 
parison is worse for the maximum infall velocity which is 



usually about 10% smaller than in the Lagrangian com- 
putation. The collapse outside of the point of maximum 
infall velocity is reproduced to an ac curacy better than 
1%. O ur protostellar results agree with lFoster fc Chevalieil 
l(l993l) to an accuracy better than 5% when we use their 
initial conditions, i.e. a 10% enhancement of pressure 
and density compared to the marginally critical Bonnor- 
Ebert sphere. We note that they do not enhance the pres- 
sure in their outer buffer region. Our simulations have 
also been compared to the slow pressure increase case of 
iHennebelle et al.l l)2003j) with a satisfactory agreement. 



2.3. Supersonic mass fraction 

Figure n shows the evolution of /s, the fraction of mass 
with supersonic infall motion, during our isothermal runs 
a nd during the simula t ion w ith the same initial conditions 
as lFoster fc Chevalieil |l993') (see section l?^ . fs depends 
critically on the initial conditions and on the parameter 
tp: the faster the increase in pressure, the more dynamical 
the collapse. This makes it a good test for the accuracy of 
the code and we reproduce the fs{t = 0) = 0.44 result of 
iFoster fc Chev"^ llT99^. 
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Fig. 1. Evolution of the supersonic mass fraction fs over 
time for different initial conditions of an isothermal col- 
lapse. Also shown is the mass fraction with motions 
greater than the surface sound speed for our reference run 
with cooling. The thin solid line is the mass over radius ra- 
tio at centre in non-dime nsional units {m/S, as defined by 
iFoster fc Chevalieilll993l) for the isothermal model with 
tp ~ 10^ yr: it is almost vertical and hence makes a pre- 
cise indicator for the instant t = 0. T he observational 
constraints set by iBelloche et al.l l(20n2j) for the class 
protostar IRAM 04191 are shown as a rectangle. 
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3. Non isothermal models with chemistry 

3.1. Numerical and physical inputs 

3.1.1. Method 

We use an improved version of a code designed to compute 
shocks in the interstellar medium (^Lesaffre et a l. 2004). 
This code makes use of a moving grid algorithm which 
resolves at best temperature and chemical gradients for 
a fixed number of 100 zones. It also solves for multi- 
fluid magneto-hydrodynamical equations fully coupled to 
a chemical network. In the present version of the code, 
we set the magnetic field to zero, wrote the equations 
in spherical coordinates and added the radiative transfer 
equations and the energy budget for the grains. 

In the following we mention only a few additional fea- 
tures that we included in the code : 

— adsorption and desorption chemistry onto and from 
grains. 

— transfer of the continuum radiation and approximate 
treatment of line cooling. 

— energy budget of the grains which transfer energy be- 
tween gas and radiation. 

Fu rther details of th e meth od can be found in lLesafFr3 
l|2002() and lAudit et"ai1 l|2002fl . 

3.1.2. Chemistry 

The chemical network used in iLesafFre et al.l l|2004j) com- 
prised 120 reactions for 32 species. This network was de- 
signed to compute the abundances of the main molecu- 
lar and atomic cooling agents of the interstellar medium. 
We added 4 species adsorbed onto grains : CO, II2O, 
CII4 and O2. We describe their ra te of adsorption onto 
grains using a model similar to iNeiad fc Wagenblasd 
1)1999(1 . Fo r the sticking coeffici e nt, w e use expres- 
sion (i) of iHollenbach fc SalDeteJ lll97j) calibrated by 
iBuch fc Zhand l)l99l|) and iMasunaga et all l|l998|) . We 
model the des orption due to cos mic-rays using the hot 
spot model of iLeger et al.l (|l984|) . The formation of H2 
molecules is the only reaction on grain surfaces included. 
Reactions of photo-dissociation, photo-ionisation as well 
as thermal sputtering of the grains are not included. Our 
treatment is thus accurate only up to grain temperatures 
of around 100 K. In the present work, we refer to the core 
as the innermost few AU where grains have such warm 
temperatures and are optically thick to radiation. We re- 
fer to the rest of the cloud as the envelope. 

3.1.3. Radiation field 

The t ransfer for the continuum radiat ion is done accord- 
ing to lDubroca fc Feugeaj l|l999|) as in lAudit et all l|2002|) 
with the Eddington factor obtained from the minimisation 
of the entropy of the pho ton gas. The Planck m ean opaci- 
ties ofthe_grains arc from lAdams fc Shil ((198(tI) (ie: similar 
to IPraine fc Le&.1983) . The outer boundary condition is 



a fixed temperature Tr — 2.7 K and the inner one is a zero 
flux. Note that by Tr we mean the holometric temperature 
of the photon gas such that the radiation energy density 
\s E = aT^ where a is the radiation density constant. 

The line transfer is assumed to be optically thin 
for the atoms and H2. This is not a good approx- 
imation deep in the core but this cooling turns out 
to be negligible ther e. The other sources of mo lecu- 
lar cooli ngs are from 'Hollcnbach fc McKe3 l(l979tl (for 
OH) and lNeufeld fc Ka,ufma,n (.19931) (CO and H2O). For 
the optical depth parame ter 7V(M) of molecule M (see 
iNeufeld fc Kaufman! I199I we use a combination of the 
column density of these molecules as well as the veloc- 
ity gradients so that radiation takes the easiest way to 
escape : 

iV(M) = n(M)/v', (1) 
with 

and 

d = R-r + A„(i?)/[1.5 X 10-2icm2 x (i?)]. (3) 

n(M) is the local density of molecule M, v is the velocity 
of the gas and Cs its sound speed, r is the distance to the 
centre, R is the outer radius, nn is the density of H nuclei 
and Ay (R) is the extinction in the visible at the outer edge 
of the cloud (A„(R) = 1 or 2 in this study ). 

As discussed bv iHoUenbach fc McKe^ l|l979|) . the es- 
cape probabilities of the photons should be modulated 
by the probability for each line that an escaping pho- 
ton be absorbed onto a grain. However, their treatment 
is only valid as long as the grains are optically thin. This 
is not valid anymore in the context of star formation af- 
ter the first adiabatic core has formed. The correct treat- 
ment would involve a multi-group treatment of the radia- 
tion field with a very fine grid in wavelength. This is still 
way out of today's computational capabilities. However, 
we tested two extreme hypotheses : either the line radia- 
tion is all absorbed by grains or the grains are transpar- 
ent to the line radiation. We found that the collapse in 
the envelope proceeds in a very similar manner for both 
hypotheses. In the following, the results are presented for 
grains transparent to the line radiation. 

Finally, the absorption of the UV field that enters only 
the expression for the heating due to the photoelectric 
effect is treated implicitly through expressions lO and © 
below. 

3.1.4. Grain temperature 

The heating and cooling processes included for the grains 
are : 

— absorption and emission of radiation by the grains (in- 
cluding the line cooling or not). The net energy rate 
gained by the grains is : 

Fradiation = aCflK{Td)nH {T^ - T^) + /3Fiinos (4) 
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where Tr, Td are the radiation and dust temperatures, 
/X is the mean weight per hydrogen nu cleus, K{Td) is 
the P lanck opacity of the grains (from lAdams 
Il986|) and c is the speed of light. Fiines is the rate of 
energy lost by the gas (hence gained by the grains) via 
emission through lines. /? = 1 when lines are assumed 
to be absorbed onto grains. [3 — Q when grains are 
transparent to the lines. 

collisional exchanges with the gas : 

rcoUisional = Ctuj^Tg^ {Tg — Td) (5) 

with a = 3.5 lO"^'* in cgs units fsee iBlackl Il987l 
page 731) and Tg is the gas temperature. 

heat released by H2 formation (2/3 or 1/3 of H2's bind- 
ing energy, depending if the line cooling is absorbed by 
the grains or not) : 



Th, =i?/Q(/3 + l)/3 



(6) 



where i?/ is the rate of formation of H2 molecules on 
grains and Q — 4.48 eV is the binding energy per H2 
molecule. 



no the grain temperature is close to the radiation temper- 
ature. However, the temperature dependence of no often 
leads to a large range of densities for which the grain tem- 
perature is close neither to Tg nor to T^ (especially for low 
Tr). For /i = 2.33 amu, k = 10"^ cm^.g-i (at Td = 10 K) 
and Td = Tg = = 10 K we get no = 6.3 10^ cm-^ 
Many authors (se e ICeccarelU eral]ll996t lGoldsmithir200lt 
iGalli et 311120021 for example) neglect rcoiiisionai in the 
grains energy budget when they estimate the grain tem- 
perature. Then they compare the radiative cooling of the 
gas to the collisional exchanges between gas and grains 
and find a critical density of 10^ cm~^ for the coupling be- 
tween gas and grain temperature, much lower than uq. It 
is hence important to consider at least all terms in equa- 
tion H10|) rather than neglecting the gas-dust collisional 
exchanges in the grains budget. 

In effect, the dust grains serve as a mediator for the 
energy transfer between gas and radiation. Indeed, plug- 
ging ((nT) into |(SJl yields the effective heating for the gas 
through gas-dust collisions expressed as a coupling be- 
tween radiation and gas temperatures: 



dust — >gas 



no 



l)-^Tg^Tr-Tg 



(13) 



— photoelectric effect on the grains is from lB1ac1<ll)l 9871) . 
page 731 : 



rphotoeioctric = 4x10 Go exp{~2.5Ay)nH 



(7) 



where Go measures the intensity of the ambient UV 
radiation field compared to the interstellar mean UV- 
field. Go is set to 1 and Ay is computed from : 



Ay{r) = Ay{R) + 1.5 X lO^^^cm^ / n^dr 

J r 



(8) 



The temperature of the grains is obtained by setting 
their energy budget to zero: 



rradiation rcoiiisionai ^¥[2 rphotoeioctric — 0. 



(9) 



Inside the collapsing cloud, the central extinction 
Av{Q) rapidly increases and shields the photoelectric ef- 
fect. If we further neglect the absorption of the lines by 
grains and the heat released by H2 formation, the energy 
budget of the grains simply reads : 

ac^i4Td)nH{T^ - T*) - anj^T^ {Td - Tg) = (10) 

The grain temperature hence lies between the gas and the 
radiation temperature. 
Equation (|10|l becomes 

no{Tr -Td) ^nniTd-Tg) (11) 

if we set 

no = ac^iKa-^T^ + iTjTr -f ST^T^ -f- T^)t; ' . (12) 

For densities much greater than no the grain temperature 
is well coupled to the gas; for densities much lower than 



3.2. Results 

As in the isothermal case, we start from a hydrostatic 
state and trigger the collapse with a slow increase in ex- 
ternal pressure. Since we now compute the temperature of 
the gas self-consistently, the parameter T of the isother- 
mal models is no longer relevant and is replaced by the 
outer boundary conditions for the radiation fields. Our 
reference run uses Tr = 2.7 K, Go = 1, Ay{R) = 1, a 
total mass of 1.7 Mq and a pressure increase time scale 
tp = 10^ yr. Another parameter is the primary ionisation 
rate by cosmic-rays C which we fixed to 5 x 10~^^ s~^. 
This parameter is important in the chem istry and for the 
heati ng of the gas through cosmic-rays l|Spitzer fc ScottI 
.1969t) . In the following, we describe the time evolution of 
this reference model, compare this model to the isothermal 
models and investigate the influence of several parameters. 



3.2.1. Description of the reference model 

Figure 13 shows the dynamical evolution of the model. 
Features such as the formation of the first adiabatic core 
and the propagation of the expansion wave where the den- 
sity profile switche s to a r~^'^ law are already described in 
the l i terature (see lLarsonlll969l: Iwhitworth fc Summers I 
ll985HMasunaga et al.lll998HMasunaga fc Inutsukall2000|) . 
We concentrate here on the heating and cooling processes. 

Figure 13 shows the evolution of the temperature pro- 
files in a sequence of snapshots during the collapse. The 
gas temperature at the very edge of the envelope is de- 
termined by the balance between the molecular cooling 
(largely dominated by CO) and the heating through the 
photoelectric effect. The temperature of the electrons is 
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decoupled from the heavier particles as they experience 
almost a ll the heating. Deeper in the envelope cosmic-ray 
heating ijSoitzer fc Scottlll969^ takes over as the UV-field 
gets shielded from the outside. 

The grains have no influence on the gas temperature 
at low densities, but the frequency of coUisional exchanges 
between the gas and the grains is proportional to the 
square of the density. When the central density rises, these 
interactions get stronger than the coUisionally saturated 
lines and CO ceases to be the main cooling agent. The 
cooling of the gas is then relayed by the grains which 
bring the gas temperature closer to the continuum radi- 
ation temperature and a dip appears in the temperature 
profile. Deeper in the core and later during the collapse the 
rate of compressional work gets stronger than the cosmic- 
ray heating and the temperature is the result of the bal- 
ance between the rate of work of the pressure and the 
coUisional cooling due to the grains. 

Finally when the grains become optically thick to the 
continuum radiation no radiative process can cool the gas 
anymore. Its temperature gets higher and the accretion 
shock takes place. We take this time as the origin of time 
for simulations with cooling. We estimate it from the in- 
stant when the maximum ofpy/p along the profile exceeds 
a half (p and are the thermal and viscous pressures) . 

Below the accretion shock the adiabatic compression 
goes on until the core temperature is high enough for the 
coUisional dissociation of H2 to proceed. The second col- 
lapse then occurs. We do not describe this phase accu- 
rately as our equation of state is for an ideal gas and we 
do not consider the evaporation of the grains. 

The code finally evolves through the singularity of the 
second collapse thanks to our fixed inner boundary. The 
density decreases in places reached by the expansion wave. 
This lowers the coupling of the gas temperature with the 
radiation temperature through the grains as seen in equa- 
tion (|13|l and the gas temperature goes up. 

3.2.2. Influence of cooling 

In this section we compare our reference run to the isother- 
mal model with same tp=10^ yr. 

We compare in figure ^ the velocity and density pro- 
files at time t = in the two simulations: cooling leads 
to much lower velocities and a less steep density profile. 
The external pressure at i = is larger for the simula- 
tion with cooling (pcxt(O) — 7.6 x 10^ K.ke.cm^'^ versus 
Pext(O) = 4.7 X 10* K.kB-cm^^), in agreement with the 
correspond i ng res ult on marginally stable states found by 
iGalli et all ^00^. Since we start both simulations with 
the same initial pext = 3 x 10** K.ke.cm"'^, the simulation 
with cooling lasts about twice as much time. 

To better understand the difference in velocities we 
trace the Lagrangian particle located at the shell of mass 
0.15 Mq in both simulations. Its velocity evolves in time 
according to the competing gravitational acceleration and 
pressure deceleration (see figure EJ. The net acceleration 
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1.00 F ' ' — ' — ' ' — ' — 1 10' 




0.01 1 , , , , , 1 10' 

100 1000 10000 
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Fig. 4. Velocity (solid lines) and density (dotted fines) 
profiles at t = for the reference simulation (thick lines) 
and an isothermal simulation with tp = 10^ yr (thin lines). 



r provides us with a typical timescale tr = \J —r jr for 
the change in radius which turns out to be of the same 
order as or lower than the velocity time scale i„ = rjr. 
We consider also ip and i^, the sound crossing time of the 
cloud (roughly = 10'* AU/0.2 km s^^ = 2.4 x 10^ yr). 
We distinguish three phases: 

1. When the pressure increase is switched on, low ampli- 
tude oscillations of period a few characterise the 
adjustment of the cloud. When they finally damp, 
the hydrostatic equilibrium is verified at all times: 

}t tp > td- Until about t = —10^ yr, the radius 
and velocity of the Lagrangian particle are hence de- 
termined by the succession of hydrostatic equilibria 
driven by the time evolution of Pcy±'- we refer to this 
phase as the quasistatic phase. 

2. For times -10*^ yr< t < -10^ yr, tp > U > The 
time scale for the variation of the radius is still greater 
than the sound crossing time and the pressure has time 
to adjust. The hydrostatic equilibrium is hence still 
rather well verified, but pext is no longer controlling 
the time scales, pext has forced the hydrostatic equi- 
librium past the marginally stable state, and the cloud 
is slowly starting to collapse. We refer to this phase as 
the detaching phase. 

3. For t > — 1 10^ yr, gets below t^. The velocity and 
radius changes are too fast for the cloud to adjust its 
pressure gradients. An increase of the gravitational ac- 
celeration refiects in a change of radius which in turn 
yields a higher acceleration. This runaway is hence sen- 
sitive to its initial conditions in radius and velocity, de- 
termined at the end of the detaching phase. We refer 
to this last phase as the collapse phase. 

We note that the transition times remain within 50% of 
these estimates for mass shells 0.3 and 0.6 M© in the case 
with cooling (with higher masses collapsing later) . But the 
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transition times can be twice earlier in the isothermal run 
compared to the run with cooling at shell mass 0.6 Mq. 
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Fig. 5. Evolution with time (-t) of the gravitational ac- 
celeration (Gm/r^ solid), the pressure deceleration (-l/p 
dp/dr, dotted) and the velocity (dashed) of the 0.15 M0 
mass shell for the reference run (thick lines) and the 
tp = 10^ yr isothermal run (thin lines). The quasistatic, 
detaching and collapse phases are indicated. 



We show a comparison of the isothermal simulation 
and the reference run in figure IHl at the very beginning of 
the collapse, ie: the last hydrostatic configuration. The col- 
lisional exchanges between gas and grains are responsible 
for a dip in the temperature profile in the inner 0.2 Mq. 
Pressure gradients depend on the density and temperature 
gradients: 



Idp ^ kB d\np ^ dT 
p dr /i dr dr 



(14) 



In the isothermal case, only the density term remains (first 
term). In the case with cooling, the temperature gradients 
(second term) are also at play. Due to the temperature 
dip, the temperature gradient term is maximum at the 
mass shell 0.05 Mq, where it is of the same order as the 
isothermal pressure gradient term, but directed outward. 
However, due to a steeper density profile in the mass range 
[0, 0.05] Mq, the density term is roughly twice the isother- 
mal pressure gradient. In effect, the resulting pressure gra- 
dient is actually of the same order as in the isothermal case 
(as noted bv iGalli ct al. 2002., for the marginally stable 
state). It is nevertheless slightly greater in the innermost 
0.05 Mq and the outermost 0.6 Mq, and slightly lower (by 
about 35% at most) in the rest of the mass. As a result, 
the bulk of the envelope (0.05 Mq < m < 1.1 Mq) is more 
extended (see figure EJ. 

This slight difference is increased during the collapse. 
Since these Lagrangian particles in the isothermal run 
start the collapse phase at a lower radius, they experi- 
ence a higher gravitational acceleration, hence they reach 



_ Temperature 

Mass 




With cooling 



4L™ 
100 



1000 

Hartiiis (AU) 



Fig. 6. Temperature of the neutral gas (solid lines) and 
mass within radius (dotted lines) at the beginning of the 
collapse (for a given density contrast of 100, around time 
t = — 9 lO'' yr) for the reference run (thick lines) and the 
isothermal simulation with tp — 10^ yr (thin lines). 



sooner lower radii and their gravitational acceleration be- 
comes even larger compared to the same Lagrangian par- 
ticle in the reference run. 

As a result the protostellar infall velocities at a given 
mass are much larger in the isothermal case than in the 
case with cooling. This explains the much lower supersonic 
mass fraction for the reference run compared to isothermal 
cases as seen in figure ^ However, to compare velocities 
at a given radius, one needs to keep in mind that a given 
mass is found at a smaller radius in the isothermal case. 
But this only partly compensates for the difference on 
the Lagrangian velocities and the net result is that infall 
velocities are larger at a given radius and a sufficiently 
late time for the isothermal run. 

Finally, we note that the reasoning "the gas cools 
down, hence there is less support, hence the collapse is 
more violent" appears to be wrong. First, the quasistatic 
sequence of equilibria approaching the marginal state is 
modified in a non trivial way by the cooling: a given exter- 
nal pressure increase yields less contraction when cooling 
is switched on. As a result, the isothermal contraction mo- 
tions are higher in the quasistatic and detaching phases. 
Second, the temperature dip in the last (unstable) hydro- 
static equilibrium leads to a more extended inner envelope. 
Third, the collapse is sensitive to the configuration (veloc- 
ities and radius) in the last hydrostatic state: it is much 
milder in the case with cooling. 

3.2.3. Influence of other parameters 

We varied tp between 10® and 10^ yr without noticing a 
significant change (the supersonic mass fraction varies at 
most by 20% in relative value). This is in contrast to the 
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isothermal models and suggests that indeed the thermal 
structure is constraining the collapse. 

We reckon that an important factor is the depth and 
steepness of the temperature dip because this regulates 
the pressure gradients in the marginally stable state. This 
is partly determined by the minimum temperature in the 
dip, which in turn depends on the external radiation field. 
We tested Gq = 0.5 which slightly lowers the temperature 
dip and yields a small decrease (of about 10% ) in the su- 
personic mass fraction. We also investigated much bigger 
changes in the value of Tr from 2.7 to 15 K while keeping 
Go set to zero. Indeed, the minimum grain temperature 
is always about 1 K above Tr when Go = 0. This pro- 
vides an easy way to control the temperature profile, even 
if the values we use do not correspond to practical situ- 
ations (remember that Tr is the bolometric temperature 
of the ambient infrared radiation field). Results show that 
velocities comparable to the infall motions in isothermal 
models are recovered for high values of T,. (above 10 K) 
which flatten the dip in temperature. 

From section 13.2.11 we note that CO is the only 
molecule that gives rise to significant cooling in the en- 
velope. In addition, the infall motions are so slow and the 
densities so high that the chemical steady state is obtained 
everywhere at all times for CO. Indeed, the adsorption 
time for the CO molecule is given by 



tad — ^dn-H 



SknTg 
nmco 



(15) 



where is the effective surface of grains per H nucleus 
and Toco is the weight of the molecule CO. On the other 
hand, the time scale for collapse can be estimated from 
the local free-fall time scale 



iff 



(16) 



where G is the gravitational constant. The ratio iff/^ad is 
greater than 1 for 



7iH > 9.4 X lO^cm-^.K/Tg 



(17) 



Since the densities at the edge of the cloud are greater than 
10^ cm^^ in the marginally stable state, the time scales for 
adsorption are significantly shorter than the compression 
time scales. This allows to compute the CO abundance 
from the equilibrium between adsorption and desorption 
onto and from grains with virtually no loss of accuracy for 
the cooling function of the gas. Since all C in gas phase is 
usually locked in CO, we conclude that out-of-equilibrium 
chemistry has no impact on the cooling in the envelope. 

Finally, we ran a simulation with a total mass of 
3.4 Mq (double of the reference run, every other param- 
eters being fixed). The density profiles are very similar 
at all times, the more massive being only an extension of 
the reference run at larger radii. As a result, the temper- 
ature dips at the beginning of the collapse are also very 
similar for radii smaller than 4000 AU (or 0.2 M©). At 
the formation of the stellar object, the innermost infall 



motions for the massive run are only about 20% higher 
than the reference model. The outermost infall motions 
are much higher (by a factor greater than 2 outside the 
shell of mass 0.2 Mq). The scaling relations expected for 
isothermal models are hence not valid for the velocity pro- 
files. 



4. Comparison with observations 

4.1. I RAM 04191 

IRAM 04191+1522 - hereafter IRAM 04191 for short - 
is one of the youngest low mass class protostars known 
so far in the Taurus molecular cloud {d = 140 pc). It 
features a prominent envelope 1.5 Mq), a powerful 
bipolar outflow and a very low bol o metric luminosity o f 
iboi 0.15 Lq ()Andre et al.lll999|) . lAndre et alJ l|l999tl 
estimated an age o f 1-3 xlO^ y r since the beginning of 
the accretion phase. iBelloche et al. (2002) performed a de- 
tailed analysis of the molecular line emission observed with 
the IRAM 30m telescope and derived strong constraints 
on the velocity structure of the envelope (see, e.g., the 
shaded area in Fig. Ej). They showed that the envelope 
is undergoing both extended infall motions and fast, dif- 
ferential rotation. They proposed that the inner envelope 
(r < 3000 — 3500 AU) corresponds to a magnetically su- 
percritical core decoupling from an environment still sup- 
ported by magnetic fields. 

Being born in the Taurus molecular cloud where star 
formation seems to proceed in a rather isolated and quiet 
manner, IRAM 04191 is likely to have experienced a slowly 
varying background in the past. It is therefore appropriate 
to compare it with models of collapse triggered quasistat- 
ically. 

4.2. Isothermal models 

Figure plots the results of two isotherma l computations 
against the observational constraints of iBelloche et alJ 
(|2002 ) for the density and velocity profiles. The two sim- 
ulations differ only by their initial conditions. 

The first simulation (figure[7^ and b) is the self-similar 
solution obtained bv .Shu (1977.) with a singular isothermal 
sphere for initial conditions. The density profiles can fit 
the observational constraints for t = to 3 10** yr but it 
shows no infall velocity in the outer part of the envelope 
for these ages. 

The second simulation (figure 0; and d) is the least 
dynamical collapse we could achieve for a Bonnor-Ebert 
sphere {tp = 10^ yr). The density profile agrees quite well 
with the observational constraints at the estimated age of 
IRAM 04191 (from lO"* to 3 x 10"* yr). But the velocities 
are much too high around a radius of 3000 AU a,tt — and 
this is even worse for the later relevant ages. The isother- 
mal collapse is therefore way too dynamical to account 
for the observed infall velocities. We plot the results for a 
temperature T = 10 K and a total mass M ~ 1.7 Mq, but 
the rescaling relations show that varying these parameters 
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Fig. 7. Comparison of isothermal models to the ob- 
servational constrai nts derived for IRAM 04191 by 
iBelloche et alJ ( 2002j) . Upper pane l fa.b): ini tial conditions 
are a singular isothermal sphere l)ShiJ ll977*>. Lower panel 
( c,d): initial conditions are a critical Bonnor-Ebert sphere. 
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Fig. 8. Com p arison of the reference model with 
iBelloche et al.l ((2002j) 's observational constraints on 
IRAM 04191. 



within their range of uncertainty (T be tween 7 and 15 K, 
M between 1 and 2 according to An dre et al.l [l9991 
has minor effects on the profiles. 

An even more stringent constraint from the observa- 
tions by Bellochc et al. (^002j is : 10% at most of the mass 
of IRAM 04191 has supersonic motions. However, even for 
the least dynamical of our isothermal simulations the su- 
personic mass fraction at time i = is still 31% and much 
more at the estimated age of IRAM 04191 (see figure 0. 

4.3. Models with cooling and chemistry 

When we include the cooling, the slower infall motions 
agree much better with the observational constraints in- 
side a radius of 3000 AU as seen in figure|Sl In particular, 
the supersonic mass fraction is now much lower and in 
much better agreement with the observations (see figure 
However, there still remains a big discrepancy outside 
this radius, where the observations require an almost con- 
stant infall velocity profile. 




t = IxlO* yr 
with cooling 





Radius (AU) 



Radius (AU) 
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Fig. 9. Com parison of mo dels with cooling at i = 10** yr 
with Bellochc et al.l ()2002|) 's observational constraints on 
IRAM 04191. 



Simulations with cooling also yield a kinetic tem- 
perature of the gas closer to the profile deduced by 
IBelloche et al.l ()2002|) from the analysis of molecular line 
spectra: in summary, they need a gas temperature of 6- 
7 K in the range of radii 2000-6000 AU and a gas tem- 
perature at the edge of around 10 K. In this respect, our 
minimum temperature (around 7 K) is too close from the 
centre of the core and our temperature at the edge of the 
core is slightly too hot (15 K). A higher Ay at the surface 
only partly remedies to the latter problem (see figure EJ: 
a bump in temperature at around 6000 AU still remains. 
This is due to the cosmic-ray heat ing. In the present work 
we use an old formulation due to ISoitzer fc Scot3 l)l969l) 
and we assume this process heats the neutral gas directly. 
However, the realistic mechanism rather involves heating 
of the electron gas which in turn exchanges heat via colli- 
sion with neutrals. This could help decrease the tempera- 
ture of the neutral gas if the coupling to electrons remains 
low. Future refinements of the code will include a proper 
treatment of the electron exchanges between the gas and 
the grains. 

The density profiles are not very sensitive to the model 
of collapse and they remain within the observational con- 
straints. Density and velocity constraints on our models 
yield an independant estimate of the age of IRAM 04191 
in the interval 0.5-2 10** yr. 



5. Discussion 

5.1. The temperature dip 

Ma,suna,ga, Tmitsukal Emh: lEva^s et a,1.l ^M); 
Zucconi et alJ l|200lll : I G alii et al.l l|2002l) also find a dip 
in the temperature profile as a result of their model for 
the grains energy budget. However, all these authors 
neglect the collisional exchanges between gas and grains 
to compute T^. As pointed out in sectional this leads to 
a coupling of and Tg as soon as the density is above 
10^ cm^'^. Besides, Td is then tied to the local radiation 
field, which is shielded near the edge of the cloud. As a 
result, the temperature starts to decrease for much larger 
radii than in our model. The temperature gradients are 
hence less steep and spread out to larger radii. Our model 
should then also account better for the flat temperature 
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profiles found outw ard of 1000 AU for the prestellar cores 
L1498 and L1517B l)Tafalla et alj r2005). 

Grain opacities and gas-grain coupling are hence the 
main physical ingredients that control the position of the 
temperature dip. Another important factor is the cosmic- 
ray heating, as discussed in section^ 

5.2. Magnetic fields and rotation 

iBelloche et al.l l|2002l) observe two regimes of infall and 
rotation in the outer envelope and inside a radius of 
3000 AU. They proposed that the outer part of the object 
is still subcritical and supported by the magnetic field, 
but the collapsing inner part is supercritical. Although we 
did not include any magnetic field, we suggest that a slow 
increase of the external pressure can mimic the slow dif- 
fusion of the magnetic field from the outer parts. Hence 
our study applies to the inner core which frees itself from 
the magnetic field. A more consistent treatment including 
a magnetic field and differential rotation is desirable but 
hard to implement in practice due to the loss of spherical 
symmetry. 

5.3. Turbulent support 

The observed macroscopic random motions amount to 
a significant fractio n of the broadening of the lines 
l)Belloche et al.ll2002|) . However, the origin of these tur- 
bulent motions is still not well understood and their role 
in the collapse is even less well known. 

5.4. Pressure variation time scale 

It is not clear whether the picture of quasistatically col- 
lapsing clouds matches the observations of the turbulent 
interstellar medium that surrounds them. Our results hold 
for pressure time scales down to tp=l Myr. However, it is 
a rather difficult task to estimate the variability of the 
turbulent background in which the cloud is embedded. 

Larson relations for a scale of 0.2 pc (corresponding 
to IRAM 04191's diameter) give a velocity dispersion of 
roughly 0.5 km.s^^. This yields a time variability of 0.4 
Myr. 

01^0(1-0) observations of iBelloche et al.l ^|2002^ with 
the IRAM 30m telescope, which probe the low density 
outskirts of IRAM 04191 (n ~ 3 x 10^ cm"^), show a 
turbulent velocity dispersion of 0.25 km.s^^ (FWHM 0.60 
km.s~^). We extrapolated this velocity dispersion from the 
beam diameter (3500 AU) to 0.2 pc with a Kolmogorov 
scaling (velocity cx lengths ) and found a time scale of the 
order of 0.35 Myr. 

These time scales are hence less than half the low- 
est value we tried for tp. However, we note that they 
are valid for the mtra-cloud velocity dispersion nowadays. 
They hence represents only a lower bound for the time- 
scale of the mter-cloud medium during the early collapse 
phase. 



6. Summary and conclusions 

We investigated the collapse of spherical clouds driven by 
a slow increase in external pressure. We compared isother- 
mal models to simulations including cooling and found 
that the sequence of hydrostatic equilibria controls the 
infall motion in both the protostellar and the prestellar 
phases. In particular, the last (unstable) hydrostatic equi- 
librium reached controls the infall velocities during the 
collapse. 

First, we note that hydrostatic equilibrium holds past 
the marginal state up to density contrasts of a factor of 10 
greater than in the marginally stable state in our reference 
simulation. A collapsing cloud still in the detaching phase 
can thus exhibit an unstable hydrostatic density profile. 
Unstable hydrostatic equilibria should then be good mod- 
els for observations of slowly contracting prestellar cores. 
Second, at the beginning of the collapse, the coupling be- 
tween the gas and the background radiation mediated by 
the grains yields a temperature dip in the inner envelope 
which slightly lowers the pressure gradients. As a result, 
the last hydrostatic configuration is more extended. Third, 
the collapse phase is very sensitive to the initial condi- 
tions, hence to the last hydrostatic state. As a result, the 
collapse is much milder than in the isothermal case. We 
therefore suggest that the external radiation field and the 
collisional coupling between gas and grains are key ingre- 
dients constraining the dynamics of the collapse. 

Our computations include a full out-of-equilibrium 
chemical network, but we show that only the adsorption 
and desorption of CO on grains have an impact on the 
cooling function. These processes can be safely considered 
at equilibrium for the range of parameters investigated. 
The energy budget is thus dominated by radiative trans- 
fer and microphysics on grains. 

We compared the results of ou r computations to ob- 
servational constraints obtained by IBelloche et alJ ( 2002?) 
for the class protostar IRAM 04191. We show that even 
the most quasistatic isothermal collapse is unable to ac- 
count for the slow infall motion observed. On the other 
hand, simulations with cooling not only improve the ki- 
netic temperature but also provide better velocity profiles 
in the inner 3000 AU. Rotation and magnetic fields are 
likely to be at play in the outer regions, but are hard to ac- 
count for in ID simulations. Future investigations should 
aim at describing ambipolar diffusion in a non-isothermal 
context. 

This study justifies the use of post-processed chem- 
istry using a prescribed dynamical background, keeping 
in mind that the line profiles directly probe infall motions 
and excitation temperatures. 
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Fig. 2. Six snapshots of the density and velocity profiles of our reference simulation with cooling. The origin of times 
is taken as the formation of the shock. 



12 



P. LesafFre et al.: The dynamical influence of cooling in the envelope of prestellar and protostellar cores. 



s 10 



Time = -2.6 lO" vr 



neutrals 
electrons 
ions 
grains 

riiilialioii 



100 
r (AU) 



Time = -7. 9 10' yr 



100 
r (AU) 



Tii:ie = -1466 yr 



100 1000 
r (AU) 



10000 



Tiinc = 53 yr 




0.01 0.10 1.00 10.00 100.00 1000.0010000.00 

r (AU) 



10000 r 



1000 r 



10 E 



Time = 198 yr 



Time = 5384a vr 




10000 F 



1000 



10 




10" 



r (AU) 



r (AU) 



Fig. 3. Evolution of the temperature profiles for the reference simulation. The times of the six snapshots correspond 
to the six snapshots displayed in figure El 
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